## resultsII.R
## Vito D'Orazio and Idean Salehyan
## script to replicate results for Who is a Terrorist
## use this script to replicate: figure 2, figure 3, figure 4, table 2, table A1, table A2, table A3, table A4, table A5

rm(list=ls())
set.seed(4422)

# to replicate results, set this to your working directory
setwd("/Users/vjdorazio/Desktop/TerrorExperiment/paper/II/Final/replication")

library(foreign)
library(readstata13)
library(Zelig) # Zelig for ordered logit
library(ZeligChoice) # Zelig for ordered logit
library(AER) # for p-values from ordered logit
library(ggplot2)

# colors
gray1<-"gray0"
gray2<-"gray50"
gray3<-"gray80"

############
## functions

f1 <- function(x1, x2) {
    return(((x2-x1)/x1)*100)
}

## functions for plotting and reworking data
oneD <- function(m, fname, model, zelig=FALSE) {
    if(zelig==TRUE) {
        mycoef <- coef(eval(m$model.call))
        myci <- confint(eval(m$model.call))
        return(data.frame(model=model, Fit=rep(fname,length(mycoef)), Variable=names(mycoef), coef=mycoef, ciL=myci[,1], ciU=myci[,2]))
    }
    mycoef <- coef(m)
    myci <- confint(m)
    return(data.frame(model=model, Fit=rep(fname,length(mycoef)), Variable=names(mycoef), coef=mycoef, ciL=myci[,1], ciU=myci[,2]))
}

cleanname <- function(d) {
    d$Variable <- as.character(d$Variable)
    d$Variable[which(d$Variable=="group")] <- "Group Only"
    d$Variable[which(d$Variable=="white")] <- "White"
    d$Variable[which(d$Variable=="white_group")] <- "White Supremacist"
    d$Variable[which(d$Variable=="arab_group")] <- "Islamic Extremist"
    d$Variable[which(d$Variable=="arab")] <- "Arab"
    d$Variable[which(d$Variable=="control")] <- "Control"
    d$Variable[which(d$Variable=="(Intercept)")] <- "Constant"
    
    return(d)
}

fill <- function(z) {
    z <- z[order(z)]
    ciL <- z[25]
    ciU <- z[975]
    mean <- mean(z)
    return(c(mean, ciL, ciU))
}

## data from the initial experiment
mydata <- read.dta13("analysisdata2.dta")

## data from the second experiment including policy recommendations
mydata2 <- read.dta("analysisdataJune20.dta")

## Table A1, visualized in Figure 2
## Models 1-2 in Appendix
fit.1 <- glm(terror ~ group + white + arab + white_group + arab_group, data=mydata, family="binomial")
fit.2 <- glm(shooting ~ group + white + arab + white_group + arab_group, data=mydata, family="binomial")


## Table A2, not visualized in the primary manuscript
## Models 3-6 in Appendix
fit.3 <- glm(shooting ~ control + white , data=mydata[which(mydata$all_group==0),], family="binomial")
fit.4 <- glm(terror ~ control + white, data=mydata[which(mydata$all_group==0),], family="binomial")
fit.5 <- glm(shooting ~ group + white_group, data=mydata[which(mydata$all_group==1),], family="binomial")
fit.6 <- glm(terror ~ group + white_group, data=mydata[which(mydata$all_group==1),], family="binomial")


## Table A3, models 9, 10, 11 visualized in Figure 3
## Models 7-11 in Appendix
fit.7 <- glm(motive_criminal ~ group + white + arab + white_group + arab_group, data=mydata, family="binomial")
fit.8 <- glm(motive_racial ~ group + white + arab + white_group + arab_group, data=mydata, family="binomial")
fit.9 <- glm(motive_religious ~ group + white + arab + white_group + arab_group, data=mydata, family="binomial")
fit.10 <- glm(motive_political ~ group + white + arab + white_group + arab_group, data=mydata, family="binomial")
fit.11 <- glm(motive_insanity ~ group + white + arab + white_group + arab_group, data=mydata, family="binomial")


## Table A4
## Models 12-13 in Appendix
fit.12 <- glm(motive_insanity ~ control + white, data=mydata[which(mydata$all_group==0),], family="binomial")
fit.13 <- glm(motive_insanity ~ group + white_group, data=mydata[which(mydata$all_group==1),], family="binomial")


## Table 2 Models 1-3 in the primary manuscript
fit.14 <- zologit$new()
fit.14$zelig(Fgun_policy ~ control + arab, data = mydata2)
tab <- coef(summary(eval(fit.14$model.call)))
p.14 <- pnorm(abs(tab[, "t value"]), lower.tail = FALSE) * 2

fit.15 <- zologit$new()
fit.15$zelig(Fhealth_policy ~ control + arab, data = mydata2)
tab <- coef(summary(eval(fit.15$model.call)))
p.15 <- pnorm(abs(tab[, "t value"]), lower.tail = FALSE) * 2

fit.16 <- zologit$new()
fit.16$zelig(Fterror_policy ~ control + arab, data = mydata2)
tab <- coef(summary(eval(fit.16$model.call)))
p.16 <- pnorm(abs(tab[, "t value"]), lower.tail = FALSE) * 2


## Table A5, visualized in Figure 4
## Models 17-19 in Appendix

fit.17 <- glm(gun_mi ~ control + arab, data=mydata2, family="binomial")
fit.18 <- glm(health_mi ~ control + arab, data=mydata2, family="binomial")
fit.19 <- glm(terror_mi ~ control + arab, data=mydata2, family="binomial")



################
#### Figure 1

plotdata <- data.frame(x=rep(c("Terrorism", "Shooting", "Unknown"),6),value=rep(0,18))
plotdata$x <- factor(plotdata$x,levels=c("Shooting","Terrorism","Unknown"))
plotdata<- plotdata[order(plotdata$x),]
plotdata$condition=rep(c("Control", "Group", "White", "White Supremacist","Arab","Islamist"),3)
plotdata$condition <- factor(plotdata$condition,levels=c("Control", "Group", "White", "White Supremacist","Arab","Islamist"))

mysums <- apply(mydata[,c("control","group","white","white_group","arab","arab_group")],2,"sum")

plotdata$value[1:6] <- apply(mydata[which(mydata$shooting==1),c("control","group","white","white_group","arab","arab_group")],2,"sum") / mysums
plotdata$value[7:12] <- apply(mydata[which(mydata$terror==1),c("control","group","white","white_group","arab","arab_group")],2,"sum") / mysums
plotdata$value[13:18] <- apply(mydata[which(mydata$shooting==0 & mydata$terror==0),c("control","group","white","white_group","arab","arab_group")],2,"sum") / mysums


d<-ggplot(plotdata, aes(y=value, x=x, fill=x)) +
    geom_bar( stat="identity") +
    facet_grid(~condition, switch="x")
d<-d+scale_fill_manual(values=c(gray1, gray2, gray3)) + theme(axis.title.x=element_blank(),legend.title=element_blank(),panel.background = element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.y=element_blank())
ggsave("figure1.pdf", d, width=14, height=7)

#####


######################
## Figures 2, 3, and 4


# some plotting arguments
pd <- position_dodge(0.5)
cols<-c(gray1, gray2)
line<-c(1,1)
mytitle<-"Model"

errorwidth<-.4
errorsize<-1.5
pointsize<-3.5

axistitle <- 19
axistext <- 16
striprel <- 1.7

temp1 <- oneD(fit.1, "1. Terrorism", "Table 1")
temp2 <- oneD(fit.2, "2. Mass Shooting", "Table 1")
out <- rbind(temp1, temp2)
out <- out[which(out$Variable %in% c("control", "group", "white", "arab", "white_group", "arab_group")),]
out<-cleanname(out)

out$Fit <- as.character(out$Fit)
out$Fit[which(out$Fit =="1. Terrorism")] <- "Terrorism"
out$Fit[which(out$Fit == "2. Mass Shooting")] <- "Shooting"

out$Variable <- factor(out$Variable, levels=c("Islamic Extremist", "White Supremacist", "Group Only", "Arab", "White","Control"))
out$Fit <- factor(out$Fit, levels=c("Shooting","Terrorism"))

mybreaks <- c("Terrorism","Shooting")

## terrorism
z5 <- zlogit$new()
z5$zelig(terror ~ control + group + white + arab + white_group + arab_group - 1, data = mydata)
z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=1, white=0, arab=0, white_group=0, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[1,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=1, arab=0, white_group=0, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[2,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=0, arab=1, white_group=0, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[3,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=0, arab=0, white_group=1, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[4,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=0, arab=0, white_group=0, arab_group=1)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[5,4:6]<-fill(fd)

## shooting
z5 <- zlogit$new()
z5$zelig(shooting ~ control + group + white + arab + white_group + arab_group - 1, data = mydata)
z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=1, white=0, arab=0, white_group=0, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[6,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=1, arab=0, white_group=0, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[7,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=0, arab=1, white_group=0, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[8,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=0, arab=0, white_group=1, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[9,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=0, arab=0, white_group=0, arab_group=1)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[10,4:6]<-fill(fd)

## ggplot for figure 2
d <- ggplot(out, aes(x=Variable, y=coef, group=Fit, colour=Fit, linetype=Fit)) + geom_errorbar(aes(ymin=ciL, ymax=ciU), width=errorwidth, position=pd, size=errorsize) + scale_linetype_manual(mytitle, values=line, breaks=mybreaks) + geom_point(position=pd, size=pointsize, colour="black") + scale_colour_manual(mytitle, values=cols, breaks=mybreaks)  + labs(x="Treatment", y="Estimated Causal Effect") + coord_flip() +  geom_hline(yintercept=0) + theme_bw()
d <- d + theme(axis.title=element_text(size=axistitle), axis.text=element_text(size=axistext), legend.title=element_text(size=axistitle), legend.text=element_text(size=axistitle), strip.text.x=element_text(size=rel(striprel)))

### saving figure 2
ggsave("figure2.pdf", d, width=12, height=7)


#############
## Figure 3
temp7 <- oneD(fit.7, "7. Crimina Motivel", "Table 4")
temp8 <- oneD(fit.8, "8. Racial Motive", "Table 4")
temp9 <- oneD(fit.9, "9. Religious Motive", "Table 4")
temp10 <- oneD(fit.10, "10. Political Motive", "Table 4")
temp11 <- oneD(fit.11, "11. Mental Illness", "Table 4")
out <- rbind(temp7, temp9, temp10, temp11, temp8)
out <- out[which(out$Variable %in% c("control", "group", "white", "arab", "white_group", "arab_group")),]
out<-cleanname(out)

out$Fit <- as.character(out$Fit)
out$Fit[which(out$Fit =="11. Mental Illness")] <- "Mental Illness"
out$Fit[which(out$Fit == "10. Political Motive")] <- "Political Motive"
out$Fit[which(out$Fit =="9. Religious Motive")] <- "Religious Motive"
out$Fit[which(out$Fit == "8. Racial Motive")] <- "Racial Motive"
out$Fit[which(out$Fit == "7. Criminal Motive")] <- "Criminal Motive"

out$Variable <- factor(out$Variable, levels=c("Islamic Extremist", "White Supremacist", "Group Only", "Arab", "White"))
out$Fit <- factor(out$Fit, levels=c("Mental Illness","Political Motive","Religious Motive","Racial Motive","Criminal Motive"))
mybreaks <- c("Crimina Motivel","Racial Motive","Religious Motive","Political Motive","Mental Illness")

out <- out[which(out$Fit %in% c("Religious Motive", "Political Motive", "Mental Illness", "Racial Motive")),]
cols <- c(gray1, gray2, gray3) # colors
mybreaks <- c("Religious Motive","Political Motive","Mental Illness", "Racial Motive")

### replicating with Zelig
## model 9 religious motive
z5 <- zlogit$new()
z5$zelig(motive_religious ~ control + group + white + arab + white_group + arab_group - 1, data = mydata)
z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=1, white=0, arab=0, white_group=0, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[1,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=1, arab=0, white_group=0, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[2,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=0, arab=1, white_group=0, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[3,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=0, arab=0, white_group=1, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[4,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=0, arab=0, white_group=0, arab_group=1)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[5,4:6]<-fill(fd)

## replicate model 10 with zelig

z5 <- zlogit$new()
z5$zelig(motive_political ~ control + group + white + arab + white_group + arab_group - 1, data = mydata)
z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=1, white=0, arab=0, white_group=0, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[6,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=1, arab=0, white_group=0, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[7,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=0, arab=1, white_group=0, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[8,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=0, arab=0, white_group=1, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[9,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=0, arab=0, white_group=0, arab_group=1)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[10,4:6]<-fill(fd)

## replicate model 11 with zelig
z5 <- zlogit$new()
z5$zelig(motive_insanity ~ control + group + white + arab + white_group + arab_group - 1, data = mydata)
z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=1, white=0, arab=0, white_group=0, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[11,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=1, arab=0, white_group=0, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[12,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=0, arab=1, white_group=0, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[13,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=0, arab=0, white_group=1, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[14,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=0, arab=0, white_group=0, arab_group=1)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[15,4:6]<-fill(fd)

## model 8 racial
z5 <- zlogit$new()
z5$zelig(motive_racial ~ control + group + white + arab + white_group + arab_group - 1, data = mydata)
z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=1, white=0, arab=0, white_group=0, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[16,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=1, arab=0, white_group=0, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[17,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=0, arab=1, white_group=0, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[18,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=0, arab=0, white_group=1, arab_group=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[19,4:6]<-fill(fd)

z5$setx(control=1, group=0, white=0, arab=0, white_group=0, arab_group=0)
z5$setx1(control=0, group=0, white=0, arab=0, white_group=0, arab_group=1)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[20,4:6]<-fill(fd)

# only plotting a subset of these
out3 <- out[which(out$Fit %in% c("Religious Motive", "Political Motive", "Mental Illness")),]
out3 <- out3[which(out3$Variable %in% c("Group Only", "White", "Arab")),]
cols <- c(gray1, gray2, gray3) # colors
mybreaks <- c("Religious Motive","Political Motive","Mental Illness")


## ggplot for figure 3
d <- ggplot(out3, aes(x=Variable, y=coef, group=Fit, colour=Fit, linetype=Fit)) + geom_errorbar(aes(ymin=ciL, ymax=ciU), width=errorwidth, position=pd, size=errorsize) + scale_linetype_manual(mytitle, values=c(1,1,1), breaks=mybreaks) + geom_point(position=pd, size=pointsize, colour="black") + scale_colour_manual(mytitle, values=cols, breaks=mybreaks)  + labs(x="Treatment", y="Estimated Causal Effect") + geom_hline(yintercept=0) + coord_flip() + theme_bw()
d <- d + theme(axis.title=element_text(size=axistitle), axis.text=element_text(size=axistext), legend.title=element_text(size=axistitle), legend.text=element_text(size=axistitle), strip.text.x=element_text(size=rel(striprel)))

## saving figure 3
ggsave("figure3.pdf", d, width=12, height=7)



##############
## Figure 4

temp17 <- oneD(fit.17, "17. Gun Control", "Table 7")
temp18 <- oneD(fit.18, "18. Healthcare for Mentally Ill", "Table 7")
temp19 <- oneD(fit.19, "19. Counter-Terrorism Policies", "Table 7")
out <- rbind(temp17, temp18, temp19)
out <- out[which(out$Variable %in% c("control", "group", "white", "arab", "white_group", "arab_group", "(Intercept)")),]
out<-cleanname(out)

out$Fit <- as.character(out$Fit)
out$Fit[which(out$Fit =="17. Gun Control")] <- "Gun Control"
out$Fit[which(out$Fit == "18. Healthcare for Mentally Ill")] <- "Healthcare for Mentally Ill"
out$Fit[which(out$Fit =="19. Counter-Terrorism Policies")] <- "Counter-Terrorism Policies"

out$Variable <- factor(out$Variable, levels=c("Constant", "Arab", "Control"))
out$Fit <- factor(out$Fit, levels=c("Counter-Terrorism Policies","Healthcare for Mentally Ill","Gun Control"))
mybreaks <- c("Gun Control","Healthcare for Mentally Ill","Counter-Terrorism Policies")
out <- out[which(out$Variable!="Constant"),]

## using zelig for simulations
## model 17 gun policy most important
z5 <- zlogit$new()
z5$zelig(gun_mi ~ control + white + arab - 1, data = mydata2)
z5$setx(control=0, white=1, arab=0)
z5$setx1(control=1, white=0, arab=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[1,4:6]<-fill(fd)

z5$setx(control=0, white=1, arab=0)
z5$setx1(control=0, white=0, arab=1)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[2,4:6]<-fill(fd)

# Model 18 health policy most important
z5 <- zlogit$new()
z5$zelig(health_mi ~ control + white + arab - 1, data = mydata2)
z5$setx(control=0, white=1, arab=0)
z5$setx1(control=1, white=0, arab=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[3,4:6]<-fill(fd)

z5$setx(control=0, white=1, arab=0)
z5$setx1(control=0, white=0, arab=1)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[4,4:6]<-fill(fd)

# Model 19 terrorism policy most important
z5 <- zlogit$new()
z5$zelig(terror_mi ~ control + white + arab - 1, data = mydata2)
z5$setx(control=0, white=1, arab=0)
z5$setx1(control=1, white=0, arab=0)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[5,4:6]<-fill(fd)

z5$setx(control=0, white=1, arab=0)
z5$setx1(control=0, white=0, arab=1)
z5$sim()
fd <- z5$get_qi(qi="fd", xvalue="x1")
out[6,4:6]<-fill(fd)

## gglot for figure 4
d <- ggplot(out, aes(x=Variable, y=coef, group=Fit, colour=Fit, linetype=Fit)) + geom_errorbar(aes(ymin=ciL, ymax=ciU), width=errorwidth, position=pd, size=errorsize) + scale_linetype_manual(mytitle, values=c(1,1,1),breaks=mybreaks) + geom_point(position=pd, size=pointsize, colour="black") + scale_colour_manual(mytitle, values=cols,breaks=mybreaks)  + labs(x="Treatment", y="Estimated Causal Effect") + geom_hline(yintercept=0) + coord_flip() + theme_bw()
d <- d + theme(axis.title=element_text(size=axistitle), axis.text=element_text(size=axistext), legend.title=element_text(size=axistitle), legend.text=element_text(size=axistitle), strip.text.x=element_text(size=rel(striprel)))

## saving figure 4
ggsave("figure4.pdf", d, width=12, height=7)

